! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2024-03-24: Initial creation.

module ffsl_common_mod

  implicit none

  character(10) :: flux_type = 'ppm'      ! Available flux types: upwind, van_leer, ppm
  character(10) :: limiter_type = 'mono'  ! Available limiter types: none, mono, pd

contains

  subroutine ppm(fm2, fm1, f, fp1, fp2, fl, df, f6)

    real, intent(in) :: fm2
    real, intent(in) :: fm1
    real, intent(in) :: f
    real, intent(in) :: fp1
    real, intent(in) :: fp2
    real, intent(out) :: fl
    real, intent(out) :: df
    real, intent(out) :: f6

    real dfl, dfr, fr

    ! Calculate values at left and right cell interfaces.
    dfl = slope(fm2, fm1, f  )
    df  = slope(fm1, f,   fp1)
    dfr = slope(f,   fp1, fp2)
    ! Why (B2) in Lin (2004) divide (dfl - df) and (df - dfr) by 3?
    fl = 0.5 * (fm1 + f) + (dfl - df) / 6.0
    fr = 0.5 * (fp1 + f) + (df - dfr) / 6.0
    ! Why (B3) and (B4) in Lin (2004) multiply df by 2?
    fl = f - sign(min(abs(df), abs(fl - f)), df)
    fr = f + sign(min(abs(df), abs(fr - f)), df)
    f6 = 6 * f - 3 * (fl + fr)
    df = fr - fl

  end subroutine ppm

  real function slope(fm1, f, fp1) result(res)

    real, intent(in) :: fm1
    real, intent(in) :: f
    real, intent(in) :: fp1

    real df, df_min, df_max

    df = (fp1 - fm1) * 0.5
    select case (limiter_type)
    case ('none')
      res = df
    case ('mono')
      df_min = 2 * (f - min(fm1, f, fp1))
      df_max = 2 * (max(fm1, f, fp1) - f)
      res = sign(min(abs(df), df_min, df_max), df)

      ! The following codes are (1.8) from Collela and Woodward (1984). It should be equivalent with the above.
      ! if ((fp1 - f) * (f - fm1) > 0) then
      !   res = sign(min(abs(df), abs(f - fm1), abs(f - fp1)), df)
      ! else
      !   res = 0.0
      ! end if
    case ('pd')
      res = sign(min(abs(df), 2 * f), df)
    end select

  end function slope

end module ffsl_common_mod